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ABSTRACT 

In this paper, by considering a simple fluid model, we investigate the role of phase 
transitions in the ISM on the galaxy scale gas dynamics. Cooling and heating timescales 
in the ISM are typically shorter than typical galactic rotation timescales, so the individual 
phases in the ISM can be assumed to be in temperature equilibrium with the radiation 
field. Using this approximation we can construct an equation of state which depends upon 
the average density and mass fractions in the individual phases. Previous studies suggest 
that there is an equilibrium phase fraction as a function of pressure. We incorporate evolu- 
tion towards this equilibrium state as a relaxation term with a time to obtain equilibrium 
t. We derive a condition in terms of a critical Mach number when one dimensional shocks 
should be continuous. For small values of the relaxation time r we show that the relax- 
ation term acts like a viscosity. We show with one dimensional simulations that increasing 
t causes shocks to become smoother. This work suggests that phase changes can strongly 
effect the gas dynamics of the ISM across spiral arms and bars. 

Subject headings: ISM: - ISM: kinematics and dynamics - ISM: shocks 



1. INTRODUCTION 

Recent studies of spiral galaxies have found that the column density of HI across spiral arms 
(e. g. in M51, Tilanus 1990), peaks at a different location when compared to the H2 (traced in 
CO) column density (Rand & Tilanus 1990). These observations, coupled with the recent work 
of Elmegreen (1993) who modeled the dependence of the molecular to atomic mass fraction in 
clouds as a function of pressure and the UV radiation field, suggest that molecular hydrogen can be 
disassociated and produced by pressure and radiation field changes across a spiral arm or galactic 
bar. In this paper, by considering a simple fluid model, we investigate the role of phase transitions 
in the ISM on the gas dynamics. Our long term goal is to construct multiphase models of the ISM 
that can be tested through observations of gas at different temperatures and densities in spiral 
and barred galaxies where changes in the pressure and radiation field are expected as a function of 
position in the galaxy. 

The interstellar medium (ISM) of the Galaxy is inhomogeneous consisting of clouds of cold 
(temperature T < 10 2 K) atomic and molecular gas, regions of warm (T ~ 10 4 K) gas, both neutral 
and ionized, and regions of hot (/ ~ 10 6 K) gas. Since the sound crossing time is smaller than 
other timescales such as the cooling timescale, it is generally assumed that the pressure in each 
phases is equal to the mean interstellar pressure, except in the molecular phase where self-gravity 
becomes important. For recent reviews see McKee (1990), Kulkarni & Heiles (1988), Scoville & 
Sanders (1987), and Begelman (1990). The first theoretical study of a multiphase ISM was the 
two phase model (Spitzer 1958, Field 1962, Field, Goldsmith & Habing 1969). Later, based on 
X-ray observations, this model was extended to include third phase by Cox & Smith (1974), and 
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McKee & Ostriker (1977) proposed that the hot phase is maintained by supernovae explosions. 
More recent studies include details of the heating and cooling balance (Shull & Woods 1985), the 
role of conductivity and cloud evaporation (Begelman & McKee 1990, McKee & Begelman 1990), 
pattern formation in conduction fronts (Elphick, Regev, & Shaviv 1992, Aharonson, Regev, & 
Shaviv 1994), and the interactions of stars and a multiphase medium in a galactic fountain (Rosen, 
Bregman, & Norman 1993, Chiang & Prendergast 1985). Most of these models have assumed that 
the average pressure in the galactic plane is static. However, the theories of spiral density waves and 
gas response in barred galaxies (e. g. Roberts, Shu & van Albada 1979) predict pressure variations 
as a function of position in a galaxy. This provides an observational setting were we can explore 
the effect of pressure changes in the ISM. 

In this paper we consider the effect of phase changes on the gas dynamics of the ISM. In 
§2 using the assumption that a equilibrium phase mix exists at a given pressure and density we 
state a simple set of continuum fluid equations which incorporate transition towards equilibrium 
as a relaxation term. The characteristics in this dynamical system are given, and a criterion in 
terms of a critical Mach number is given for when shocks should be continuous. In §3 a set of 
one dimensional simulations are presented. It is shown from these simulations that increasing the 
relaxation time causes the shocks to become more smooth. For a small relaxation time, r, we show 
that the relaxation term acts like a viscosity. A short discussion follows in §4. 

2. Continuum Equations 

In §2.1 we show that cooling timescales are generally faster than rotational timescales in a 
galaxy. This justifies using a first approximation where each phase is considered to be in thermal 
equilibrium with the radiation field. This approximation is particularly desirable since dynamical 
equations to describe the thermal state of the fluid in each phase are not required (see Pistinner 
& Shaviv 1993 for a more general set of multiphase fluid equations). We can then use continuum 
fluid equations where the variables represent averages of local quantities over size scales larger than 
individual clouds, and far fewer equations are required to describe the multiphase fluid. In §2.2 
we describe a simple set of equations to describe such a system. We model the evolution of the 
phase mix in terms of evolution towards an equilibrium state with r, the relaxation timescale or 
timescale required to obtain equilibrium. In §2.3 in the one- dimensional system sound speeds are 
given for two limits of the system, r — ► and r — ► oo. To gain insight on this set of equations in 
§2.4 we derive a dispersion relation for perturbations from equilibrium. We show that for small r, 
the relaxation term causes damping similar to that caused by viscosity. 

2.1 Timescales 

We estimate the cooling timescales in order to compare these timescales to galactic timescales. 
In the cold gas T ~ 100°K, Pottasch, Wesselius & van Duinen (1979) from observations of the CII 
158/i line estimate that the cooling rate per Hydrogen nucleon is 1.2 X 10~ 25 erg/s. This implies 
a cooling timescale of t coo i ~ kT/(dE/dt)~ 10 4 years which is fast compared to galactic scales. 
Likewise, in the ionized medium with T ~ 8000° K and the cooling for the same pressure is about 
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the same (Kulkarni k Heiles 1988, Shull k Woods 1985) so that kT/(dE/dt)~ 4 X 10 5 years. For 
comparison, one galactic rotation period in a galaxy with circular rotation speed v c at a radius r is 

2vrr 7 /200km s~ 1 \ f r \ , 

r = _„, xl „ v (2 . L1) 

Variations in pressure and density caused by a spiral density wave or barred potential should extend 
over a solid angle which corresponds to a time that is larger than few percent of the rotation 
period. Therefore, except in the inner regions of galaxies (e. g. in nuclear rings, where the rotation 
timescale is faster) heating and cooling timescales are generally faster than the gas response to a 
barred potential or a spiral density wave. This implies that in each phase, the state of the gas can 
be described as a function of the radiation field and the average pressure. We note that a warm 
(T ~ 8000° K) neutral component is estimated to have a much longer cooling time ~ 10 7 years (see 
references in Kulkarni k Heiles 1988) and so may not be in equilibrium with the radiation field. In 
this paper we do not consider the role of this component of the ISM. 



2.2 Continuum Fluid Equations 

Using the approximation that each phase is in thermal equilibrium with the radiation field, we 
do not require equations to describe the thermal state of the fluid in each phase. We adopt here a 
formalism where the variables are averages over size scales larger than the size of individual clouds. 
We assume here that the mean velocity v of all phases is the same so that all phases move together. 
Mass conservation is then given by 

^ + V-(pv) = (2.2.1) 

where p is the average density summed over a scale greater than that of individual clouds. Conser- 
vation of momentum or the Euler's equation is 

<9v VP 

_ + (v . V )v=-- (2,2.2, 

where P is the pressure, which is a function of the phase mix as well as p. 

Previous studies suggest that at a given pressure there exists an equilibrium phase mix. Begel- 
man k McKee (1990) derived an equilibrium equation of state as a function of pressure by con- 
sidering the role of conduction for size scales larger than the Field length on the phase mix, and 
Elmegreen (1993) considered the equilibrium molecular to atomic gas fraction as a function of the 
pressure and radiation field. If we assume that the medium requires a relaxation timescale r to 
achieve equilibrium, which we expect is greater than the cooling timescale in the individual phases, 
then we can describe the evolution towards equilibrium as a relaxation process. For /, the mass 
fraction of one phase of a two phase medium, we can describe the evolution towards equilibrium 
as 

^kll + V • (pfv) = 9 { jeqm ~ f) (2.2.3) 

where r is the timescale required to obtain equilibrium and we have written the equilibrium phase 
fraction as f eqm . We expect f eqm to be a function of local variables such as the pressure and 
the radiation field. We note that relaxation terms can strongly effect the dynamics of the Euler 
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equation and can produce dissipation like terms in this equation (Whitham 1974 in traffic flow and 
molecular gases, Glimm 1986 in multiphase fluids, Liu 1991, Landau & Lifshitz 1959 in molecular 
gases, see bulk viscosity). 



2.3 Sound Speeds and Critical Mach Number 

We note that our set of equations is similar to the equations that describe relaxation effects 
in gases, where the internal energy may lag behind the equilibrium value. For example vibrational 
and rotational degrees of freedom may take longer to reach equilibrium than kinetic or translational 
ones in a molecular gas (see references cited above). In this section we follow the intuitive discussion 
by Whitham (1974) to derive a criterion for when shocks are continuous rather than sharp followed 
by a region of relaxation. 

We can view the set of previous equations in terms of two systems. If we consider timescales 
of order r then we require the full set of equations listed above; we call this system shock system I 
(Whitham's language). The characteristics for this system are v , v ± c for v the velocity and c the 
sound speed where 

c 2 = dP< g> f \ (2.3.1) 

If we average the fluid flow over timescales larger than r then we can consider the phase fraction 
to be at the equilibrium value. This is equivalent to the limit r — ► and is described by only two 
equations (neglecting the equation for the phase fraction), however the phase fraction is then equal 
at all times to the equilibrium value and the characteristics are v ± c', with sound speed c' where 

c' 2 = c 2 - d 2 (2.3.2) 

and 

dP df eqm QP 

A 2 _ 9/ dP d P , 2 _ 3 _3^ 



dP df 

In general we expect 4^ > 0. For / representing the denser phase fraction we expect that at higher 

pressure there should be a higher fraction of dense gas so that d ^ 9 p m > 0. Similarly we expect 
< when the averaged density, p, is kept constant. Therefore for reasonable functions d 2 > 

and 

v-c<v-c'<v<v + c'<v + c (2.3.4) 

which is a general condition for stability, in other words if the fluid begins with a smooth initial 
condition, discontinuities will not develop. 

Viewed from the system II (averaging over timescales larger than r) which takes the flow 
between two uniform states in equilibrium, the shock will be continuous in system I if 

v 2 - c 2 < U < mi + ci (2.3.5) 

for U the velocity of the shock. Subscripts here refer to quantities before and after the shock. A 
frozen shock followed by a region of relaxation will occur when 

U>u 1 +c 1 (2.3.6) 
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or when the Mach number M 

U — u-i c . 
M = > -. 2.3.7 

c\ c 

The critical Mach number M cr n = c/c'. 



2.4 Dispersion Relation, and Dissipative Behavior 

Sometimes a dispersion relation can be used to gain some insight into the stability and dissipa- 
tional behavior of a set of dynamical equations. Here we derive the dispersion relation for our set 
of equations and show for small r that the flow is stable and dispersive. We derive the dispersion 
relation by considering the first order perturbations p = po + ep a , v = vq + ev a , and / = / + e/ a 
where p a ,v a , f a ,<x e lkx ~ lu ' t an d e is small. We obtain the following dispersion relation 

(-^ + 2kv u> + e (c 2 + vl)) (-iu + tkv + i - i^p) + y^r^^ = (2- 4 - 1 ) 

for c given above. In the limit r — ► 0, u = kv^ ± A;c' and in the limit r — ► oo, a; = fcfo ± A;c. The 
solutions are stable and dispersive when all roots of the dispersion equation have Im(uj) < 0. For 
small t we can expand the dispersion relation in powers of r and to first order in r we find that 

ik 2 T 8P df eam 8P 

u = ^ ±fa + __Jg=._. (2 .4.2, 

As described in the previous section, we we expect ^ > 0, d ^jT > and < so that the 
imaginary component of a; should be less than zero, and the solutions are damped. If we incorporate 
a kinematic viscosity v into equation (2.2.2) and neglect phase changes, then we find to first order 
in v that 

u = kv ± kc - (2.4.3) 

It is clear from comparison of these two previous equations that the relaxation term in equation 
(2.2.3) behaves for small r similar to a viscosity. If the ratio of the dissipation terms 



dp df eqm op 

df dP dp 



> 1 (2.4.4) 



then the dominant form of dissipation is due to phase changes for a system that has both phase 
changes and viscosity. This inequality should be roughly equivalent (in order of magnitude) to 



rc 2 

— > 1 (2.4.5) 



which we expect is satisfied when r is greater than the mean collision time of clouds. 
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3. One-Dimensional Simulations 

In this section we describe the results of a set of one- dimensional simulations of the equations 
given in §2.2. We run a one- dimensional computer code that is a 3rd order (in Li norm) finite 
difference scheme in conservative form which uses an ENO-LLF (Essential Non-Oscillatory, Local 
Lax-Friedrichs) shock capturing scheme (Shu & Osher 1988, 1989). The time stepping method 
is a 3rd order TVD (total variation bounded) Runge Kutta method. For a description of the 
implementation see Don & Quillen (1995) and Quillen (1995). Only the left hand side of equation 
(2.2.3) can be put in conservative form. We therefore treat the right hand side of this equation as a 
forcing term. Equations (2.2.2) and (2.2.1) are easily put in conservative form. Our simulations are 
run with 128 grid points with periodic boundary conditions in a two phase model. As described in 
the appendix we assume that the pressure is proportional to a power of the density in each phase 
(see equation A. 5 for the form and notation). The denser phase is described with parameters A\, 
pi and ai, whereas the less dense phase is described with parameters A2, P2 an( i a 2- We use a 
model where a\ = a.2 = a. We assume an equilibrium phase fraction function (e. g. Elmegreen 
1993) 

fegm(P) = B (3.1) 

but don't allow f eqm to be greater than 1. Parameters for the simulations are listed in Table 1 and 
are unitless. The simulations run with grid spacing dx = 1 and are updated at a time interval of 
dt = 0.2. 

We run a series of simulations with different values of r and initial Mach number. Initial 
conditions are listed in Table 2 for the individual simulations which are denoted S a , Sb, T a and T5. 
Each simulation begins with a shock velocity that is zero in our inertial frame, and a phase fraction 
on either side of the shock that is at the equilibrium value. Initial conditions are a step function 
with initial values pi, fi and vi one the left side of the interval and p r , f r and v r one the right side 
of the interval. Because of the periodic boundary conditions the shock at the edge of the interval 
is a rarefr action wave. 

In Figure 1 we present two simulations (denoted S a and Sb) with Mach number M below the 
critical Mach number (see §2.3) that have different relaxation times r. S a has r = 0.5 and Sb has 
t = 5.5. Note that the simulation with larger r (Sb) has smoother shocks. This is because of the 
dissipative effect of the relaxation term in equation (2.2.3). 



4. SUMMARY AND DISCUSSION 

By incorporating phase changes in a simplistic way into a continuum fluid equation, we find 
that phase changes can significantly effect the gas dynamics of the ISM. The sound speed for the 
system with phase equilibrium should be lower than the sound speed of a system with a frozen 
phase mix. For small values of r the timescale required to obtain an equilibrium phase mix, phase 
changes act similar to a viscosity. Simulated one- dimensional shocks with larger values of r are 
smoother than those with smaller values of r, consistent with the dissipative nature of the system. 

Many physical processes have not been considered in this paper. The role of radiation which 
can affect the phase mix and possibly cause star formation mediated radiative precursors may also 
strongly affect the continuum fluid dynamics. Turbulence and magnetic fields, also not considered 
here, could dominate the pressure (in the Euler's equation). In the centers of galaxies where rotation 
timescales are faster, (and in the outer parts of galaxies at small scales) departure from thermal 
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equilibrium may occur. At small scales the effects of individual clouds cause the continuum fluid 
approximation to be invalid. It has been suggested (e. g. Quillen et al. 1995) that different phases 
in the ISM may not be moving together. In future more complicated sets of dynamical equations 
which can incorporate these processes should be studied. 

It has recently become possible to observe the ISM in galaxies in a variety of wavelengths 
(and so phases); for example, it may become possible to constrain the radiation field and the 
cooling rate of the cooler components of the ISM as a function of position in a galaxy using far 
infrared observations. We hope that multi-wavelength observations of galaxies are coupled with 2- 
dimensional (and eventually 3-dimensional) fluid simulations of these galaxies, which can be run in 
a variety of fluid systems. This type of study will make it possible to constrain the fluid properties 
of the ISM. 

We note that in barred and spiral galaxies the process of dissipation may determine the shape 
and location of the shocks as well as the gas inflow rate (Roberts, Shu & van Albada 1979, Barnes & 
Hernquist 1991). Whereas galactic disks are rather inefficient accretion disks (Steiman-Cameron & 
Durisen 1988) when shocks are formed in a non-axisymmetric gravitational potential mass transfer 
can occur rapidly (e. g. Mihos & Hernquist 1994a, b, Shlosman & Noguchi 1993). One possible 
explanation for this is that although the sheer viscosity may be small, the bulk viscosity (which 
becomes important in shocks and is present in numerical simulations) may be large. We note 
that when relaxation effects are important, it is primarily the bulk viscosity that is increased 
(e. g. Landau k Lifshitz 1959). In the barred galaxy NGC 7479, Quillen et al. (1995) have 
measured a gas inflow rate along the bar. It is likely that gas inflow in barred galaxies is strongly 
dependent upon the nature of the fluid (ISM). We hope in future to explore how gas inflow rates 
may depend upon the constituents of the ISM. 

We acknowledge helpful discussions and correspondence with C. Thompson, B. Elmegreen, R. 
Oiling, G. Shaviv, 0. Regev, N. Shaviv, J. Glimm, A. Gould, K. Griest, P. Schecter, A. Goodman 
& D. DePoy. A.C.Q. acknowledges the support of a Columbus fellowship. 
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Appendix. Equation of State, an Example 

We can write the average density p summed over a scale greater than the scale of individual 
clouds as 

P = Y,fv i Pi (A.l) 



P 1 = Y^fmiPi 1 

i 

fm, P 



fv, 



Pi 



(A.2) 



(A.3) 



where f m . , /„. and pi are the mass fraction, volume filling fraction and density of phase i respec- 
tively. We can also write 

/> = J> (A.4) 

i 

for G{ = f m> p the mass per unit volume of phase i averaged over a volume which contains all phases. 
For an opticaly thin tracer of a particular gas phase the intensity observed would be proportional 
to G{ integrated along the line of sight. 

The two phase model assumes that the pressure in all phases are the same. Studies of heating 
and cooling balance find that the pressure P and njj the number density of H atoms are related 
in each phase by scaling laws (Shull & Woods 1985) that look like 



P 



gas 



(A.5) 



for various assumptions about the UV, X-ray and cosmic ray fluxes. Using this we can write the 
density in each phase as 

Pi = n H p t = (P/Ai)^ a 'p t (A.6) 

for a mean molecular weight in phase i of pi. These equations provide enough information to solve 
for an "equation of state", P(p, f m ,)- This is is easiest to do if a is the same for all phases, in which 
case 

-1 A/a 



P( P J m ,) = p C 



(A.7) 



The derivatives ^ and gj^— are needed to do simulations and estimate sound speeds, and these 
derivatives can be written 



dP 
dP 

dfm, 



E 



fv, 
(Xi 



P 



Pi 



E 



fvj 



(A.8) 



(A.9) 



We note that the molecular gas is normally not included in the "three phase model of the 
ISM" (McKee 1990) because most molecular clouds are self-gravitating (Larson 1981, see McKee's 
review) and so are not in pressure equilibrium with the rest of the ISM. However, if there is a 
molecular to atomic phase transition, the pressure is indirectly dependent upon the mass fraction 
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in the molecular phase. As Begelman & McKee (1990) have noted, physical processes causing phase 
changes such as conductivity and cloud evaporation can significantly effect the mean pressure and 
energy density of the medium. We can incorporate molecular clouds into the above model by 
considering their integrated density as a function of the pressure at the surface of the clouds. 
Elmegreen (1993) by considering statistics of self gravitating and diffuse molecular clouds, and self 
shielding of these clouds in a background UV radiation field, finds that the equilibrium molecular 
to atomic gas fraction is proportional to a power of the pressure. 
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TABLE 1 
Parameters for One-D Simulations 



Parameter 



A 1 


1.0 


A 2 


2.0 


Mi 


1.0 


Pi 


1.0 


a 


1.3 


B 


0.4 


(1 


1.0 



TABLE 2 

Initial Conditions 
for One-D Simulations 1 



Parameter S a ,Sb T a ,Tb 

M l7l3 

M crit 1.31 

Pi 1-32 

p r 0.86 

vi 0.85 

v r 1.29 

fi 0.80 

/ r 0.60 



1 M is the initial Mach number and M cr a is the critical Mach 
number (see equation 2.3.7). pi, vi and /; are the initial 
density, velocity and mass fraction of the dense phase on the 
left hand side of the shock. p r , v r and f r are the initial 
density, velocity and mass fraction of the dense phase on the 
right hand side of the shock. 
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FIGURE CAPTIONS 



Figure 1. Comparison of two M < M cr a shocks with different values of r, the relaxation time. 
Parameters are listed in Tables 1 and 2. The solid line is simulation S a and the dotted line is 
simulation Sb- The two simulations have identical initial conditions except that S a has relaxation 
time t = 0.5 and Sb has relaxation time r = 5.5. The simulations are shown at a time of / = 11.8. 
Note that the one with larger r has a smoother profile, consistent with the dissipative behavior 
caused by the phase changes. 



